Edge properties and Majorana fermions in the proposed chiral rf-wave superconducting state of 

doped graphene 
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We investigate the effect of edges on the intrinsic d-wave superconducting state in graphene doped close to 
the van Hove singularity. While the bulk is in a chiral d x 2_ y 2 + id xy state, the order parameter at any edge is 
enhanced and has d x n_ y i -symmetry, with a decay length strongly increasing with weakening superconductivity. 
No graphene edge is pair breaking for the d x 2_ y 2 state and there are no localized zero-energy edge states. We 
find two chiral edge modes which carry a spontaneous, but not quantized, quasiparticle current related to the 
zero-energy momentum. Moreover, for realistic values of the Rashba spin-orbit coupling, a Majorana fermion 
appears at the edge when tuning a Zeeman field. 

PACS numbers: 74.20.Rp, 74.70.Wz, 74.20.Mn, 73.20.At, 71.10.Pm 



Graphene, a single layer of carbon, has generated immense 
interest ever since its experimental discoverylllj . Lately, ex- 
perimental advances in doping methods J2, 0] have allowed 
the electron density to approach the van Hove singularities 
(VHSs) at 25% hole or electron doping. The logarithmically 
diverging density of states (DOS) at the VHS can allow non- 
trivial ordered ground-states to emerge due to strongly en- 
hanced effects of interactions. Very recently, both perturbative 
renormalization group (RG) and functional RG calcula- 
tions |2j,|6[] have shown that a chiral spin-singlet d x 2_ y 2 +id xy 
(d±+id,2) superconducting state likely emerges from electron- 
electron interactions in graphene doped to the vicinity of the 
VHS. This is in agreement with earlier studies of stron g int er- 
actions on the honeycomb lattice near half-filling lASS- 

Rather unique to the honeycomb lattice is the degeneracy of 
the two d-wave pairing channels lu ll 10 . Below the supercon- 
ducting transition temperature (T c ), this degeneracy results 
in the time-reversal symmetry breaking d\ +id 2 state [1, 4]. 
However, any imperfections, and most notably edges, might 
destroy this degeneracy and generate a local superconducting 
state different from that in the bulk. At the same time, many 
of the exotic features prop osed for a di +id 2 superconductor, 
such as spontaneous 112, T^, or even quantized ffl. edge 
currents and quantized spin- and thermal Hall effects 
are intimately linked to its edge states. In order to determine 
the properties of d±+id 2 superconducting graphene, it is there- 
fore imperative to understand the effect of edges on the super- 
conducting state. 

In this Letter we establish the edge properties of di+id 2 su- 
perconducting graphene doped to the vicinity of the VHS. We 
show that, while the bulk is in a d% + id 2 state, any edge will 
be in a pure, and enhanced, di-wave state. Due to a very long 
decay length of the edge di state, the edges influence even 
the properties of macroscopic graphene samples. We find two 
well-localized chiral edge modes which carry a spontaneous, 
but not quantized, edge current. Furthermore, we show that by 
including a realistic Rashba spin-orbit coupling, graphene can 
be tuned, using a Zeeman field, to host a Majorana fermion at 
the edge. These results establish the exotic properties of the 
chiral di+id 2 superconducting state in doped graphene, which 



if experimentally realized, will provide an exemplary play- 
ground for topological superconductivity. Furthermore, these 
results are also very important for any experimental scheme 
aimed at detecting the d\ + id 2 state in graphene, as such 
scheme will likely be based on the distinctive properties of 
the edge. 

We approximate the 7r-band structure of graphene as: 



Hn = -t 



E 



c icr C i(7) 
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where t — 2.5 eV is the nearest neighbor (NN) hopping 
amplitude and Ci„ is the annihilation operator on site i with 
spin a. The chemical potential is \i and the VHS appears at 
|U = ±t, where the Fermi surface transitions from being cen- 
tered around K, K' to T. We study two different models for 
superconducting pairing from repulsive electron-electron in- 
teractions: 
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In the limit of very strong on-site Coulomb repulsion (mean- 
field) pairing appears on NN bonds such that a a = 5 a 
(a = 1, 2, 3) [ID, whereas a moderate on-site repulsion gives 
rise to pairing on next-nearest-neighbor (NNN) bonds with 
o Q = 7 a J2l, see Fig. [TJ a )- The high electron density near 
the VHS efficiently screen long-range electron-electron in- 
teractions, and we also show that our results are largely in- 
dependent on the choice of a. In mean-field theory the or- 
der parameter can be calculated from the condition A a (i) = 
-J{ciiCi +aa t - c^c i+ a a i)- Here J is the effective (con- 
stant) pairing potential arising from the electron-electron in- 
teractions and residing on NN bonds for a = 5 and on NNN 
bonds for a — 7. Using this condition for A, the Hamilto- 
nian H = Hq + can be solved self-consistently within 
the Bogoliubov-de Gennes formalism J3,[l8t]. The favored 
bulk solution of A Q belongs to the two-dimensional E 2 irre- 
ducible representation of the Cq v lattice point group. This rep- 
resentation can be expressed in the basis a c i 1 = (1, — h, — h), 
which has d\ = d x 2_ y 2 symmetry when Hq is diagonal, and 
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hd 2 = (0, 2>~ ^r) w hich has d 2 = d xy symmetry, see 
Fig. [Tib)- In the translational invariant bulk, these two so- 
lutions have the same T c , but below T c the complex combina- 

I — I r-i 

tion di+id.2 has the lowest free energy QJJ, |4J] . There is also an 
s-wave solution, a s = (1,1,1), but it only appears subdomi- 
nantly and at very strong pairing. 
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FIG. 1: (Color online) (a) Graphene with NN bonds 5 a , NNN bonds 
y a , zigzag and armchair edges indicated, (b) Different d-wave super- 
conducting order parameters for NN and NNN pairing with negative 
(blue) and positive (red) sign. 

In order to quantify the edge effects we study thick ribbons 
with both zigzag and armchair edges. We assume smooth 
edges and Fourier transform in the direction parallel to the 
edge. Due to computational limitations we need J > O.bt 
in order to reach bulk conditions inside the slab. This gives 
rather large A Q , but by studying the J-dependence we can 
nonetheless draw conclusions for the experimentally relevant 
low- J regime. 

Superconducting state at the edge. — In the bulk, the d\+id2 
state has a free energy AF lower than the dx,2 states, which 
varies strongly with both doping and pairing potential, see in- 
set in Fig. [2c). However, sample edges break the translational 
invariance and a qualitatively different solution emerges. Fig- 
ure |2ja) shows how the zigzag edge completely suppresses 
the imaginary part of A a , while at the same time enhancing 
the magnitude. This suppression leads to a pure di solution 
at the edge, an effect we quantify in Fig. |2jb) by plotting the 
dx-character (Ai - ±(A 2 + A 3 ))| 2 . The edge behavior 
can be understood by noting that bonds <5 2 and ^3 (72 and 73) 
are equivalent for both armchair and zigzag edges IU9II and, 
therefore, the di-wave state is heavily favored at both type of 
edges. Since the edge is of the zigzag type for edges with 
30° and 90° angles off the x-axis and of the armchair type 
for 0° and 60° angles, we conclude that any edge should host 
d-wave order with nodes angled 45° from the edge direction. 
In order to quantify the spatial extent of this edge effect, we 
calculate a decay length £ by fitting the d\ -character profile 
to the functional form (C'e~ x/ ^ + 0.5) with C « 0.5. As 
seen in Fig. Etc), £ varies strongly with AF, but very little 
with edge type and doping level. Furthermore, the increase 
in £ for NNN pairing compared to NN pairing suggests that 
the edge will be even more important in models with longer 
ranged Coulomb repulsion. The strongly increasing £ with 
decreasing AF has far-reaching consequences for graphene. 
For example, J — 0.5t and doping at the VHS gives £ ss 25 A 
for NN pairing. With an expected much weaker superconduct- 
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FIG. 2: (Color online) (a) Order parameter profile for the zigzag edge 
for NN J = 0.75t at the VHS with real (black) and imaginary (red) 
part for Ai (thick), A2 (thin), and A3 (dashed) [black dashed line 
is hidden behind black solid line since Re(A2) ~ Re(A3)]. (b) 
Character of the order parameter in (a): d\ (thick black), d^ (thick 
red), and s (thin black). Dotted line marks the bulk value, (c) Decay 
length £ of the di -character as function of AF for different doping 
levels, edges, and superconducting pairing: NN pairing, zigzag edge, 
and /i = t (black x), fi = 0.8t (red o), /1 = 1.2t (green o) or 
armchair edge and /1 = t (blue x ), NNN pairing, zigzag edge, and 
fi = t (black *), fj, = O.St (red □), fi = 1.2* (green A) [blue, 
x symbols is often completely overlaying black, x symbols since 
no notable difference is found between zigzag and armchair edges]. 
Inset shows AF as function of the pairing potential for NN pairing 
(black) and NNN pairing (red) for fx = t (thick) and fj, = 1.2t (thin). 
fi < t has a AF curve similar to fi > t. 



ing pairing in real graphene, the edge will not only modify the 
properties of the superconducting state in graphene nanorib- 
bons, but also in macroscopically sized graphene samples. We 
have verified that both the d±+id2 state itself and edge effects 
described here are stable in the presence of random disorder 

iH. 

Chiral edge states. — Any di + ic?2 state, even with one sub- 
dominant part, violates both time-reversal and parity sym- 
metry and has been shown to host two chiral edge states 
H2I ffl [l5ll . The topological invariant guaranteeing the ex- 
istence of these two chiral edge modes also causes quantized 
spin- and thermal-Hall responses Ell. Figure |3 a) shows 
the band structure for a zigzag slab. The self-consistent solu- 
tion (thick black) gives two Dirac cones located at ±fco, where 
bands with same velocities reside on the same surface, thus 
yielding two co-propagating chiral surface states per edge. 
The band structure for the constant (non-self-consistent) bulk 
di+id2 state also has two Dirac cones (thin black), but shifted 
away from fcrj. The shift is directly related to the d\ state at 
the edge. The d\ state has no surface states on the zigzag 
edge, only bulk nodal quasiparticles, where the nodes for a d\ 
order parameter with amplitude equal to that on the edge are 
located at ±fco (thin red). The similarity between the d\ +id2 
and g?i edge band structures thus makes for only modest ef- 
fects of the edge on the self-consistent band structure. It also 
results in the chiral edge modes being well localized to the 
edge, as seen in the local density of states (LDOS) plot in 
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FIG. 3: (Color online) (a) Band structure for a zigzag edge slab with 
NN J = 0.75t, fi — t, and self-consistent A (thick black), constant 
di+id,2 state corresponding to the bulk state (thin black), and constant 
di state corresponding in amplitude to the di state at the surface, (b) 
LDOS across the ribbon for the self-consistent solution in (a) inter- 
polating between 0.2 (black) to (white) states/eV/unit cell, showing 
a bulk gap of 0.18 eV and gapless edge states, (c) Quasiparticle edge 
current in units of e/h as function of superconducting bulk order pa- 
rameter A(l, e 27 " /3 , e ^ i/z ) for zigzag edge with ^ = t (black x), 
/i = 0.8t (red o), and armchair edge with \i = t (green A). 



Fig. Ob). The constant edge LDOS is a consequence of the 
one-dimensional Dirac spectrum. We note especially that no 
d-wave superconducting graphene edge will display a zero- 
bias conductance peak due to zero-energy surface states, in 
contrast to the cuprate superconductors [1311 . Such a peak is 
only present when the order parameter for incidence angle 9 
on the edge has a different sign from when the angle is tt — 9. 
This only happens for the d,2 -solution on both the zigzag and 
armchair edge. 

The breaking of time-reversal symmetry gives rise to spon- 
taneous edge currents carried by the chiral edge modes 11121 - 
By combining the charge continuity equation with 
the Heisenberg equation for the particle density [4], we calcu- 
late in Fig. [3|c) the quasiparticle edge current as function of 



H x = i\R ^2 z- (s CT!<7 / x dij^Cj^ , (3) 

where dy is the unit vector from site j to i. Supercon- 
ducting two-dimensional systems with Rashba spin-orbit cou- 
pling and magnetic field have recently attracted much atten- 
tion due to the possibility of creating Majorana fermions at 
vortex cores or edges J5, 2j] 22]. At edges the Majorana 
fermion appears as a single mode crossing the bulk gap. This 
should be contrasted with the behavior found above, where 
the edge instead hosts two modes. We will here show that a 
Majorana mode is created in d-wave superconducting doped 
graphene in the presence a moderate Zeeman field: Hh = 
-h z J2i(4t Ci t ~ C 4 C 4)- Due to spin-mixing in H\, the ba- 
sis vector = (c^c^Ci-fC^) has to be used when express- 
ing the Hamiltonian H cx t = Hq + Ha + H\ + Hh in ma- 
trix form: H ext = ^X^T~L cxt X . This results in a doubling 
of the number of eigenstates compared to the physical band 
structure. This doubling is necessary for the appearance of 
the Majorana fermion, since a regular fermion consists of two 
Majorana fermions. 

A change in the number of edge modes marks a topolog- 
ical phase transition which, in general, can only occur when 
the bulk energy gap closes. We therefore start by identify the 
conditions for bulk zero energy solutions of H cxt . Close to the 
VHS we can, to a first approximation, use only the partially 
occupied 7r-band for small A, A#, and h z . A straightforward 
calculation [5] for this one-band Hamiltonian gives the fol- 
lowing bulk-gap closing conditions at ^ ~ t: 



(v-t\e k \f 



|A fe ||Afl£fc| =0, 



(4) 



where ek = ^2 a e zkSc " is the band structure, ipk = arg(efe), 
Afc = — ^ Q A a cos(fc<5 a — (fk) is the fc-dependent intra- 
band superconducting order for NN pairing [1], and Ck = 
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of the bulk order parameter A(l, e 
evidence for a quantized boundary current equal to 2eA//i, as 
previously suggested [14J. In fact, we find a non-linear rela- 
tionship between current and A, a strong variation with dop- 
ing level, and, most importantly, the armchair current even 
decreases when A increases. The last result can be under- 
stood by studying the zero-energy crossing ±fco of the chiral 
edge modes. For the zigzag edge fco increases with increasing 
A, whereas for the armchair edge fc decreases. In general, 
we find that changes in current are proportional to Sky with 
P « 1 — 2. This, at least, partially agree with earlier results 
reporting a /3 — 2 dependence 11211 . Finite fc-point sampling 
and neglecting the screening supercurrents could potentially 
explain the discrepancy. 

Majorana mode. — Heavy doping of graphene, by either ad- 

I — L ^"L 

atom deposition [2] or gating [3], breaks the z — ¥ —z mirror 



4 W3) We find no ^ s t ^ le spin-orbit interaction when expressed in the form H\ = 



Cka' for the one-band model. Equa- 



tions (01 are met at T, K, and M in the Brillouin zone, where 
they produce the conditions (/i — 2>t) 2 — h 2 z , fj 2 — h 2 ., and 
(fi - tj 2 + Afc(M) = h 2 z , respectively. At [i — t only 
the last condition is satisfied for small h z , which is neces- 
sary for superconductivity to survive. We find Afc(M) = 2 A 
for the A(l, e 27 ™/ 3 , e 47 ™/ 3 ) order parameter and, thus, at the 
VHS there is a topological phase transition at h c = 2A. Fig- 
ure Ufa) shows how the eigenvalue spectrum of a supercon- 
ducting zigzag slab at the VHS develops when h z is swept 
past h c . At finite A/? and/or h z the chiral modes in Fig.|3la) 
split with one mode moving towards k y = and the other 
one towards the zone boundary at k y = tt, see left-most figure 
in Fig. Sta). At h c (center figure) the bulk gap closes at both 
k y = 0, tt. The closure at k y = tt annihilates the outer chiral 
modes whereas the closure at k v = leaves a new Dirac cone 
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crossing the bulk band gap with the two modes belonging to 
different edges. Thus, at h z > h c we are left with three modes 
per edge crossing the bulk gap. The odd number establishes 
the existence of a Majorana mode alongside the two remnant 
chiral modes. Figure HJb) shows how A develops in the pres- 
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ing state in heavily doped graphene is in a pure d\ state on 
any edge. The d\ edge state significantly modifies the super- 
conducting state even in macroscopic graphene samples due 
to a long decay length. Moreover, d\ + ids superconducting 
graphene hosts two well-localized chiral edge modes, which 
carry a non-quantized spontaneous quasiparticle current. A 
Majorana mode can also be created at the edge by tuning a 
moderate Zeeman field. These results establish the properties 
of the d\+id% state in graphene, and will be important for any 
experimental detection of this state. 

The author thanks A. V. Balatsky, M. Fogelstrom, and 
T H. Hansson for discussions and the Swedish research coun- 
cil (VR) for support. 



FIG. 4: (Color online) (a) Eigenvalue spectrum for a zigzag slab with 
NN J = 1.2*, n = t, X R = 0.2*, and h z = 0.4, 0.535 and 0.6 eV 
(left to right), with h c — 0.535 eV. Small gaps in the surface states 
are due to limited fc-point sampling, (b) Self-consistent A as function 
of h z for J = 1.2* (black), 0.9* (red) for X R = 0.05* (thick), 0.2* 
(thin), and 0.3* (dashed). Dotted line mark the h c — 2A one-band 
model phase transition. Crosses mark the numerical phase transition, 
(c) Eigenvalue amplitude squared for the Majorana mode in (a) for 
h z = 0.54 eV (black), 0.56 (dashed), and 0.6 (red). 

ence of an applied Zeeman field h z , with x -symbols marking 
the phase transition into the phase with a Majorana fermion. 
The dotted line marks the one-band result h c = 2A, which is 
a good approximation for small Xr. In this small A/j-regime 
there is a very pronounced drop in A at the phase transition 
with only a small remnant superconducting state in the Ma- 
jorana phase at h z > h c , which results in a poorly resolved 
Majorana mode. Larger A^ gives a stronger superconduct- 
ing state in the Majorana phase. However, for Xr > 0.2* we 
find h c > 2A, and the superconducting state is again very 
weak beyond the phase transition. We thus conclude that, in 
order to create a Majorana fermion at the edge of <i-wave su- 
perconducting graphene doped very close to the VHS, a small 
to moderate Rashba spin-orbit coupling, Xr ~ 0.2*, and a 
Zeeman field of the order of 2A is needed. With reported 



tunability with electric field |24], as well as impurity-induced 
Rashba spin-orbit coupling 12511 . Xr ~ 0.2* is likely within 
experimental reach in heavily doped graphene. The Zeeman 
field can be generated by proximity to a ferromagnetic insula- 
tor, whereas if applying an external magnetic field, orbital ef- 
fects also needs to be taken into account. Finally, in Fig.|4fc) 
we plot the spatial profile of the Majorana mode amplitude 
just beyond h c . Due to the larger A at the edge, the bulk 
enters the Majorana-supporting topological phase before the 
edge. Therefore, the Majorana mode does not appear at the 
edge but is spread throughout the sample for h z > h c . Not 
until h z > 2 A (edge) does the Majorana mode appear as a 
pure edge excitation. 

In summary, we have shown that the di+id2 superconduct- 
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Supplementary material 

In this supplementary material we provide: (1) a detailed, largely self-contained, description of the method underlying our 
results, and (2) numerical data showing the relative robustness of the bulk d\ +id,2 superconducting state and its edge properties 
in the presence of disorder. 

Method. — As described in the main text, we use the Hamiltonian H = Hq + Ha, where 

H = -t c L c j> + mX c ^ Cict ' (5) 
h a = J2 A ^)(4 t 4 + a ai - 4A +aat ) + H.c. (6) 

i,ct 

Here Cj CT is the annihilation operator on site i of the honeycomb lattice with spin a, t = 2.5 eV is the nearest-neighbor (NN) 
hopping amplitude, and fj, is the chemical potential, where fj, = ±£ corresponds to the van Hove singularities (VHSs) (H is 
particle-hole symmetric so hole and electron doping give the same result). Furthermore, the superconducting order parameter 
A Q resides on NN bonds when a a = 5 a and on next-nearest neighbor (NNN) bonds when a a = "f a , where a — 1,2,3 labels the 
three inequivalent bond directions, see Fig. 1(a) in the main text. Within mean-field theory, A a is defined by the self-consistent 
condition: 



A„(i) = — J{cn,Ci + a a -f - c^c i+aa i), (7) 

where J is the effective pairing potential on NN bonds for a = 5 and on NNN bonds for a = 7. J is a consequence of the local 
repulsive Coulomb interaction, which in the limit of very strong on-site repulsion results in pairing on NN bonds [ 1 ] whereas a 
moderate on-site repulsion results in NNN bond pairing Q. 

We can solve H within the Bogoliubov-de Gennes formalism by writing 

H = X^UX with X^ = (c[ t ,c a ), (8) 

and diagonalizing the matrix H to find all eigenvalues E v and eigenvectors V ', where v = 1, 2N for N sites. We can then 
define new operators = (7J) using X = VY with the columns of V given by the eigenvectors V v , such that the Hamiltonian 
H is diagonal in these new operators: H = E v rf v ^ v . A self-consistent solution scheme start with first guessing the value 
of A Q , diagonalizing H for this value, using the self-consistent condition Eq. (O to recalculate A Q from the eigenvalues and 
eigenvectors, and then reiterate these steps until the order parameter A Q within two subsequent steps changes less than a small 
predetermined convergence limit. Using the self-consistent value for A Q any electronic property of the system can be calculated 
in using the eigenbasis. For example, the local density of states (LDOS) can be calculated as 

Di(E) = X \Vf\ 2 S(E - E v ) + \V£ +i \ 2 S(E + E»), (9) 

where the first part is the spin-up contribution and the second part the spin-down contribution. Numerically, we use a small 
Gaussian broadening for the 5-functions. We are also interested in the quasiparticle current, which can be calculated using the 
continuity equation for the charge current density J: 

V-J + |j=0 (10) 
at 

together with the Heisenberg equation for the particle density per unit cell nc 

drii i , 

where p — e(n) QSl- The quantum average of the commutator in Eq. (fTTT > can easily be shown to only contain Hq for a 
self-consistent solution of A Q . The total quasiparticle edge current is then simply / = J11, where the summation is over all 
unit cells at the edge with a finite J parallel to the edge. 

The above formalism can be applied to any structure on the honeycomb lattice. In order to investigate edge properties, we 
study H on thick graphene ribbons having either zigzag or armchair edges. We make sure that the ribbons are always thick 
enough to guarantee bulk conditions in the interior. For simplicity, we assume smooth edges so we can Fourier transform in the 
direction along the edge, which introduces a fc-point index, while reducing the site index i to only enumerate sites perpendicular 
to the edge, i.e. % now only measures the distance to the edge. 
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We also study the influence of a finite Rashba spin-orbit coupling: 

H x = i\n ^2 z- (s CT!<7 / x dy)cj CT c,>' , (12) 

(i,j),a,cr' 

where is the unit vector from site j to i, in combination with a Zeeman field: 

H h = -h z ^( c lt c *T " c i4. c 4)- ( 13 > 

i 

The spin-mixing in the Rashba term now requires us to write H cxt = Ho + Ha + H\ + Hh as 

tfcxt = ^H cxt X with Xt = (c^atcu), (14) 

i.e. we need to double the number of eigenstates compared to the physical band structure. This is expected since a Majorana 
fermion is essentially half a fermion. By applying the same self-consistent procedure described above to H cxt , we can solve for 
the superconducting order parameter A, calculate all physical observables such as the LDOS, and also calculate the eigenvalue 
spectrum, which contains the Majorana mode for a large enough field h z . The Majorana mode appears when the eigenvalue 
spectrum develops from having an even number of edge modes to an odd number. Such a change in the number of edge modes 
is in general always associated with the closing of the bulk gap. We can analytically extract the approximate bulk gap closing 
condition from an effective one-band model. The kinetic Hamiltonian Hq is diagonalized in the bulk by changing the basis from 
the site-operators {ca, cb} on the two inequivalent sites A and B, to the band operators {a, b} through: 

( C A ka \ _ J_ ( a ka + b ka \ 

I, c Bka ) J2\ e-*V" (a ka - b ka ) ) " KLD) 
Here a} k(y creates an electron in the lower 7r-band and creates an electron in the upper 7r-band, such that 



H = ^ [(P _ t£k)a\ a a ka + (ji + te k )b\ a b ko 

ka 



(16) 



where the fc-dependence of the 7r-bands is given by = | ^ Q e lk ' Sa | and ip k — arg e. tk Sa )- We will for simplicity now 
assume fi ~ t and focus on the lower 7r-band, but the same calculation is also valid for doping levels around the VHS at /i = —t. 
By only keeping terms within the lower band, ignoring effects in the upper band along with any cross-terms, we arrive at: 



a L ak ° 



H' z = -h z a 

H' A = - ^Aa cos(/c • 5 a - ¥>k){alt a -ki) + Hc - 

k,CX 



kacr' 



where L k = A i? Im[ e -^ fc (-^ e tt * 2 + ^e^ 3 , e lkSl - \e lkS2 - \e lkS3 , 0)]. The one-band Bogoliubov-de Gennes Hamiltonian 
H' ext = H'o + H' A + H' z + H' x can now be diagonalized and we find the eigenvalues 



E{k) = y (M - ^fe) 2 + + K + |A fc p ± 2^J(fi - te k f\%Ll + [(/x - ie fc ) 2 + \A k \*]hl (18) 

where A k = J2a cos(fc • S a — tp k ). Following the procedure in Ref. J[ the zero energy values of Eq. dl81 l. or equivalently 
the bulk-gap closing condition, satisfy 

(fi-t\e k \) 2 + A 2 k = h 2 + X 2 R \£ k \ 2 , \A k \\\ R £ k \ = : (19) 

which is the same as Eq. (4) in the main text. From this equation we locate the only low-field bulk closing point to be at 
(ji - t) 2 + A\{M) = h 2 , where A k (M) = 2A for the A(l, e 27Ti/3 , e 47 ™/ 3 ) order parameter. As seen in Fig. 4(b) in the main 
text this approximative bulk closing condition is accurate for small to moderately large \r. 
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Disorder effects. — Heavy doping of graphene will undoubtedly introduce some amount of disorder into the system. Disorder 
can affect the results derived in this Letter in several ways. First of all, sufficiently strong disorder will suppress the super- 
conducting order parameter, this is especially true in non- s-wave superconductors. Secondly, disorder breaks the translational 
invariance and, thus, the two d-wave channels in graphene are no longer guaranteed to be degenerate. Related to this is the fact 
that there exists also an extended s-wave solution, which, in general, is heavily disfavored but in the presence of disorder might 
become more important. In addition to these bulk effects, disorder might also influence the edge properties of the di+id,2 state. 

In order to study the effect of disorder we model both the bulk and zigzag edges in the presence of Anderson disorder, i.e. a 
locally fluctuating chemical potential: 

#o,dis = -t ^2 Ao c i<y + + <^;) C L C ^ ( 2 °) 

where the local chemical potential variations 6 fa are distributed randomly within the interval [— W, W), with W being the 
disorder strength. This type of disorder model captures the effect of local charge inhomogeneities introduced by the doping. It 
is also reasonable to assume, as is done here, that the disorder will in general be directionally independent, such that it does 
not single out one bond direction over the other. We solve i/o,dis + Ha for multiple disorder configurations in large bulk and 
edge samples and study how the superconducting order evolves with the disorder strength W . We fix fi — 1 which maximizes 
the effect of the disorder, since both negative and positive deviations from the VHS causes the superconducting order parameter 
to decrease. In Fig. [2c) we plot on the right axis the superconducting order parameter Ai as function of the NN pairing 




FIG. 5: (Color online) (a) Order parameter Ai profiles at the zigzag edge for NN J = 1.2t at the VHS for a 20 unit cell wide sample with 
W = 0.5t = 1 eV disorder (red). Clean sample (thick black), (b) Average character of the order parameter in (a): di (thick black), dz (thick 
red), and s (thin black). Dashed lines display the clean di and d-2 results, respectively, (c) Bulk Ai order parameter (right, black axis) as 
function of NN J for clean (black dotted), W = 0.2t — 0.5 eV (black solid), and W = 0.5i = 1 eV (black dashed) and the corresponding 
s-wave character (left, red axis). 

potential J for a 40 x 40 A bulk sample. The results are averaged over as many as 20 different disorder configurations. For 
W = 0.2i = 0.5 eV there is a suppression of the superconducting state for J < 0.4t, at which point the character of the 
superconducting state also changes from perfect d\ +id% to contain a significant amount of s-wave character (right axis). At 
J = OAt, W is 18 times larger than Ai and thus the di+id2 state survives disorder at least an order of a magnitude stronger than 
the superconducting gap. The appearance of a sizable s-wave component at very strong disorder is expected since isotropic states 
are more robust against disorder, but this also suppresses the overall superconducting order parameter. For W = OAt = 1 eV 
we find the same scenario, with the d\ +id2 state being suppressed into a weaker partial s-wave state at J < 0.6t, where W 
is 12 times larger than Ai. Based on these results, we expect the d\ + id2 state to survive essentially unchanged in the bulk in 
the presence of even moderately strong disorder. In Figs. |3a,b) we plot the behavior at the edge for a representative W = 1 eV 
disorder configuration. Computational demands limit the size of the sample and we are forced to use a rather large J = 1.2t. 
Nonetheless, the disorder strength is still in this case almost 3 times larger than the bulk Ai value. The Ai profile into the sample 
in Fig. [51 a) shows a noticeable spatial variation, but still, the average is not suppressed much from the clean limit. As seen in 
Fig.[5jb), the average character of the superconducting state is also essentially left unchanged by this relatively strong disorder. 
We thus conclude that even moderately strong disorder does not influence the edge properties of the d\ +id2 superconducting 
state in heavily doped graphene. 



[1] A. M. Black-Schaffer and S. Doniach, Phys. Rev. B 75, 134512 (2007). 

[2] M. Kiesel, C. Piatt, W. Hanke, D. A. Abanin, and R. Thomale, arXiv: 1109.2953 (unpublished). 

[3] L. Covaci and F. Marsiglio, Phys. Rev. B 73, 014503 (2006). 

[4] A. M. Black-Schaffer and S. Doniach, Phys. Rev. B 78, 024504 (2008). 

[5] M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B 82, 134521 (2010). 



